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ABSTRACT 


The trailing vortices generated by the control planes of sub- 
marines give rise to surface signatures in the form of scars and 


Striations. 


Two counter-fou eae vortices were generated in a novel experi- 
mental system and their interaction with the free surface was investi- 
gated. In addition, the governing equations have been solved through 
the use of the boundary-element method for a representative Froude 
number. The results have been expressed in terms of the depth of 
submergence of the vortices, their mutual induction velocity, and the 
initial vortex separation. It has been shown that the free surface 
begins to deform when the vortices are at a distance of about one ini- 
tial vortex separation from the free surface. The height of the maxi- 
mum deformation is attained at a normalized time of about 0.1, when 
the vortices are at a distance of about 0.5 bo from the free surface. 
The elevated part of the surface is bounded by two scars, whose 


motion is slaved to that of the vortices. 
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I. INTRODUCTION 


Vortices and vortex wakes have become a major theme of 
aerodynamics research since the advent of the large aircraft and the 
understanding of their evolution required an examination of many of 
the fundamental problems in fluid mechanics. Much of the progress 
made during the past two decades was discussed at the Symposium on 
Aircraft Wake Turbulence and Its Detection [Ref. 1] and the Aircraft 
Wake Vortices Conference [Ref. 2]. Comprehensive reviews of the 
entire subject have been given by Donaldson and Bilanin [Ref. 3], 
Widnall [Ref. 4], and Hallock and Eberle [Ref. 5]. 

These studies, as well as numerous others carried out since 1977, 
have uncovered a number of complex problems which must be 
resolved in order to achieve a better understanding of the important 
features of trailing vortices in homogeneous and stratified media. The 
principal ones are as follows (Ref. 6]: 


1. Roll-up process: The velocity and turbulence distribution at any 
station behind the wing depend on the wing section, wing-tip 
shape, Reynolds number, wing incidence, and the distance of 
the station from the wing [Ref. 7]. The distributions of the initial 
velocity and turbulence, which influence the roll-up and the 
decay process, cannot be changed independently. For example, 
a change in the tip shape changes the core size, as well as the 
velocity and turbulence distributions. High levels of turbulence 
result in an increased diffusion of vorticity, which in turn 
increase the core size. 


2. Probe sensitivity of the vortices: Flow visualization studies 


suggest that trailing vortices are extremely sensitive to 
disturbances created by even very small probes or bubbles. This 
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forces one to use non-intrusive means of meaSurement such as a 
Laser Doppler Velocimeter. Even then, “vortex wandering” 
[Ref. 8], which makes the vortices appear larger than normal in 
time-averaged velocity measurements (for vortices generated by 
a wing in a wind tunnel), or the unsteady nature of the flow (for 
vortices generated by a wing in a tow basin) makes the mean 
velocity profiles in the vortices difficult to determine. 


Large-scale instabilities: The vortices are seldom observed to 
decay away owing to viscous and turbulent dissipation, but are 
almost always destroyed by either mutual induction instability 
(Crow instability [Ref. 9]) and/or vortex breakdown. The Crow 
instability grows exponentially, and results either in linking of 
the vortex pair into a series of crude vortex rings or in a highly 
disorganized intermingling of the vortices. 


Vortex breakdown, whose mathematical details have not yet 
been adequately treated, rearranges the vortex structure and 
increases the core size, turbulence, and energy dissipation. 
Thus, it is very difficult to measure accurately the trajectories of 
the three-dimensional vortices from their creation to their 
ultimate demise. 


Reynolds number: Even the highest Reynolds numbers, based . 
on wing chord, reached in wind tunnels or towing basins, are an 
order of magnitude lower than what is possible for an aircraft. 
Thus, the scale effects are not easy to assess. 


Ambient conditions such as turbulence and stratification play 
major roles in the evolution of vortices. The quantification of 
these effects requires numerical analysis and extremely careful 
experiments [Ref. 10]. 


Ground or free surface effects: The vortex pair may move toward 
a rigid boundary at which the no-slip condition must be satisfied 
or toward a free surface at which the zero-shear condition must 
be satisfied. In either case, the vortices come under the 
influence of their images and move accordingly. 


The phenomenon is further complicated by several additional 


facts. When the vortices are propelled toward a rigid surface, vorticity 


of opposite sign is generated on the no-slip boundary and swept 


toward the vortex pair. The total vorticity diminishes very quickly as 


vorticity from the two regions diffuses, the wall region serving as a 
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strong sink for the vorticity associated with the original vortex 
[Ref. 11]. The development of a boundary layer along the rigid wall 
may give rise to flow separation for sufficiently high Reynolds numbers. 
With or without such a separation, however, the center of the vortex 
pair eventually moves away, or “rebounds,” from the wall [Refs. 1l- 
13]. a = 

For the case of a zero-stress boundary, the free surface still acts as 
a vorticity sink, but this is relatively weak due to the absence of 
intense oppositely signed vorticity. Thus, in the absence of other 
impending phenomena, one expects a mild interaction between the 
vortices and the free surface and a small rebound of the vortex pair 
from the free surface. However, the ability of the free surface to 
deform under the influence of strain fields leads to a strong inter- 
Betion between the vortices and the free surface. 

It is evident from the foregoing that the motion and the life-span 
of trailing vortices are governed by a number of nonlinearly dependent 
complex phenomena. A number a experimental and analytical studies 
have been carried out at the Naval Postgraduate School by Sarpkaya 
and his students [Refs. 6, 14-21] in order to investigate the effects of 
‘these parameters on the rise and demise of the trailing vortices in 
homogeneous and density stratified media. These studies have clearly 
identified the various demise mechanisms in both media and estab- 
lished basic relationships between the rise of vortices and the 
governing parameters in a finite as well as effectively infinite medium 


[Ref. 20], free from ambient turbulence. 


15 


The present investigation is a continuation and refinement of the 
previous studies. The intent is to analyze the rise and demise of a 
vortex pair in a medium with a deformable free surface. This problem 
is of interest both from the standpoint of determining the interaction 
effect of the free surface on the rise and demise of the vortex pair and 


from the standpoint of predicting the resulting free surface shape- 
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II. PHYSICAL ANALYSIS 


A. DIMENSIONAL ANALYSIS 
The dependent parameter of major importance to the problem 
solution-is the instantaneous position of the vortex pair (x, y). It may 


be expressed as a function of the following parameters [Ref. 22]: 


meertit, Vo, do, Oo, AP/dy, V, Do, g. Te, &, L113) le) 
and 
y = f(t, Vo, do, Po, dp/dy, v, Do, g, re, €, L11) (2) 


in which the variable definitions are as follows: 


I! time 

Vo initial mutual induction velocity of the vortices 

dg initial depth of the vortex pair 

Do reference density of the medium 

dpo/dy linear density gradient 

Vv kinematic viscosity of the medium 

Do initial separation of the vortex pair 

g gravitational acceleration 

Te effective core radius of the vortex 

E rate of decay of the turbulent energy per unit mass 


Lii integral scale of the turbulent field 
The height and width of the test section were not included in the 
foregoing because a detailed analysis, based on ideal vortices, has 


shown that the velocities induced by the bottom or sides were 
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negligible. Effects of surface tension on the instantaneous position of 
the vortex pair are deemed negligible and thus are not included as a 
parameter in equations (1) and (2). 


A dimensional analysis of equations (1) and (2) yields: 
X/Do = f(Vot/Do, do/Do, NoBo/Vo. Vo2/gDbo, Vobo/V, Te/Do, €*, Ly 1/Do) (3) 
and 


y/Do = f(Vot/bDo, do/Bo, Nobo/Vo. Vo2/gbo, Vobo/V, re/Do, €*, Li1/Do) (4) 


in which 
_ Y2 
jes ts S| (5) 


is known as the Brunt-Vaisala Frequency. 

All parameters in equations (3) and (4) may be changed indepen- 
dently except re/Do, which is taken as nature provides it. The primary 
reason for this is that a century of theoretical and experimental aero- 
dynamics research has been incapable of describing the details of the 
structure of the tip vortex to be used as the initial conditions in the 
viscous solution. It is surprising, but true, that until recently the 
importance of the generating surface shape and its influence upon 
both the initial tangential velocity profile and the initial turbulence in 
the vortex core had not been fully appreciated. Here the said influ- 
ence has been characterized in terms of an effective core radius with 


full awareness of its shortcomings. 
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B. GOVERNING DIFFERENTIAL EQUATIONS AND BOUNDARY 
CONDITIONS 


The generation of internal waves and the rise and demise of a 
vortex pair in a stratified medium may be analyzed through the use of 
the equations of motion for an incompressible fluid. These equations 
may be applied for both laminar and turbulent motions provided that a 
suitable turbulence closure model is used and the usual Boussinesq 
approximation (gravitational acceleration is much larger than fluid 
acceleration) is adopted. For the type of motion considered herein, 
the Boussinesq approximation is quite valid and has been used in the 
investigation of all types of internal waves in stratified fluids. 

For a two-dimensional flow, with y vertical and x horizontal, the 


Navier-Stokes equations of motion are: 


dus dues—é—* 1 0P 
att UOx + Voy = —p oxt vV7u (6) 
OV dV ov 1 OP 9 


Differentiating equations (6) and (7) with respect to y and x respec- 





tively yields 
du d2u_ss—s o2uu—s—s dvou—s—sén2tuu_—si1 Op OP's ~02P OF es 
dyat * Oxdy eUineen a, + By oy + Ya2y = pZ oy ax p oxoy * Yay V u) = (8) 


ov dudv dv dv dv _ 1 dp oP 
dxot ' oxox “Ox2 * Oxoy * “axdy ~ p2 ox oy 





Subtracting equation (9) from equation (8) yields 


Subtracting equation (9) from equation (8) yields 


2@-)-30G-2]-306-8]-232- 


1 0p OP of out _ Ov 

o2 ox ay t YY €: ox (10) 
Seles 6) for + & yield 7 
olving equation (6) for 9 ox vie S 

1 oP D 

5 TW U- De (11) 
where 

D() oa), a), oa) 

Dt = t + Ux TV Oy 


Solving equation (7) for ~ $ yields 


D 
yet WV - Br (12) 


Substituting the results of equations (11) and (12) into equation (10) 
and defining vorticity as 


C= (55 oe) 13) 


20 


results in the following expression: 


AB on $B Go BR) 


=, (14) 


— 


When the gravitational acceleration is several orders of magnitude 
larger than the fluid accelerations, the Boussinesq approximation is 
appropriately invoked and the terms in brackets in equation (14) may 
be neglected. 

In order to deal with the case of a density stratified medium, the 
density is defined as 


P =Po+ Ply) + p'(x.y,t) ~ (15) 


where fo is the reference density, B(y) is the initial variation of density 
with y, and p'(x,y,t) is the fluctuating part of the density with time. 
Then 


dp _ dp’ 
ix ox 


and equation (14), neglecting the bracketed terms, becomes 


26 Bul), BW yay £20 a6 
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The last term above represents the effect of the density gradient and 
gives rise to oppositely signed vorticity in a nonhomogeneous fluid. 


The diffusion of density is given by 


P+ ugh + ve = vV%p (17) 


Substituting equation (15) for density in the above equation yields — 

20" uO. o(B 2) = (SEE SE ag) 
The equation of continuity is 

ou av 


2 =o (19) 


Adding equation (19) to the left side of equation (18) and simplifying 


results in 
- 5 
Be tae ae a Vay + Ie ve 7 


Equations (16) and (20) are thus the governing differential equa- 
tions for the motion of a vortex pair in a density stratified fluid. 

For the development of boundary conditions, it is assumed at the 
outset that the influence of the sidewalls and bottom of the test 
section on the rise and demise of the vortex pair is negligible. The 
validity of this assumption, when used in conjunction with numerical 
computations, has been demonstrated by previous investigators [Refs. 


18-20]. The fluid domain can thus be considered to be bounded only 
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by a free surface which can be described as a function of x and time as 


follows: 
y = (x,t) 
where 
N(x,0) = O 


describes the initial undisturbed location of the free surface. 

The free surface requires both a kinematic and a dynamic bound- 
ary condition as described by Sarpkaya and Issacson [Ref. 22]. The 
kinematic condition states that any particle which lies on the free 


surface at any instant will never leave it. This leads to 
Dn _ on, one . 
fm ct Vox “VY at y=n (21) 
The_dynamic free surface condition requires that the pressure 
difference across the interface results in a force normal to the bound- 
ary which is due wholly to surface tension. This condition takes the 
form 


1 l 
P =P, + ofr * Ra) (22) 


1 1 
where o is the surface tension, 55) and Ro are the radii of curvature of 


the free surface in any two orthogonal directions, and P - Pg is the 


pressure difference across the interface. 
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As noted previously, the effects of surface tension can be assumed 
to be negligible. This observation has been verified experimentally by 
Gray [Ref. 17]. Thus, if the flow is considered inviscid, the pressure P 
within the fluid can be described by the unsteady Bernoulli equation as 


follows: 

@ Ls gn +o = Fit (23) 
where © is the velocity potential such that 

u=s— Vv =5p (24) 


and q? = u2 + v2. F(t) is an arbitrary function of time only. 
When the pressure just outside the liquid is constant (i.e., atmo- 


spheric), the free surface condition reduces to 


eB) 


Za) 2 
e+e en=0 at yan (25) 


where F(t) has been included as part of eo 
C. DIMENSIONLESS PARAMETERS 

It is convenient at this point to cast the governing equations, (15) 
and (20), in nondimensional forms, scaling each variable by a quantity 
characteristic of its expected magnitude. Two possible time scales 
exist. The dynamic time scale utilizes the time a characteristic length 
would be traversed by a fluid particle traveling at the characteristic 


velocity. The buoyant time scale is based on the natural buoyancy 


24 


frequency of the stratified flow, i.e., the Brunt-Vaisala frequency No 
defined by equation (5). Each of these scales gives a slightly different 


form of the normalized governing equations. 


1. Dynamic Scale 


Introducing U, and Le as the characteristic velocity and 


length, one has the following nondimensionalized quantities: = 
Cm = CL¢/Ue tm = Uct/Le Um = u/Ue Vm = v/Uc 


Xm = X/Le Ym = y/Le Pm = P/Po (26) 


Substituting the above into equation (16) yields 


Uc _d (GmUc), 1 9 U —- — chs U CmUc) _ 
Tedtm\ Le f* le dxm|"™~e Tg Gi, (WEES AT |S 








rz va (=) + ose (Pop'n) ee 
Simplifying 

Ud dm , fs Basa) mis 2a) PS V2 tn + +i em 
which becomes 
or 


where 


Re = UcLc/v (characteristic Reynolds number) 
Fy = Us/vV glc 


(characteristic Froude number) 


Similarly, the density diffusion equation may be expressed in 


nondimensional terms by making equivalent substitutions into equa- 
tion (20) as follows: 


' - Pde 
Ps eae) oye) = -v wap! + vse 


This becomes 


3 1 Qa l a 

Te Stim (PoP ia) + Lo Sig UctimPoPin) * Te By WmUeP opal 
‘i a2 

7 .— (Pop'm) + yz m (PoP 'm) + [Z i 


= Yin 


Which upon simplification reduces to 


UcPo 0'm | UcPo UMP) ri UcPo O(VmP'm) 
Seu Oth) Slaeno eee 








om 
_UcPo  OPm = ev, VPo 070m 
“Le “™ Gym LZ Ym Pm + 72 Gye 


Multiplying both sides by Le/poUe yields 
0P'm z (Um?) | o(VmP'm) 


dp Pm | Vv 9 2 070m 
dtm Om Om  °™ ym  Uele [va Pmt Oy2- | (30) 
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summaries 


Incorporating the characteristic Reynolds number, this equation 


becomes 
00'm OuUmP'm  AVvmPm) — APm Ili, , 070m 
dtm? oxm * dym ~ “™ dym* Re|’mPmt Gz (31) 


The Brunt-Vaisala frequency given by equation (5) as 


_|=& \* 
No = [58 2 


may be written as 


N72 = -g 0(Po Pm) 
07 Pole dm 


which becomes 


Noble am 


g ~ dm 





To further simplify equation (31), the following definitions are made: 


N2 = g/Le (characteristic Brunt-Vaisala frequency 
squared) (33) 


and 





= 2. 22m | (34) 


Substituting equation (34) into equation (31), the nondimensional 


form of the density diffusion equation in final form becomes 
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(35) 


m 


OP'm 9UmPm) | AVmPm) _ lficjo. ,?Pm 
im *  oxm * Oym 7 “m+Re[VmPm+ oz 


2. Buoyant Scale 
Again introducing U,; and Le as the characteristic velocity and 


length, the following nondimensional quantities are defined: 


u Vv Xx 
Um = Uo VEL > Wie cit ia ym * Lo Pm *5 


In this case, however, time is nondimensionalized using the square of 


the characteristic Brunt-Vaisala frequency. Namely 


NZ = g/Le 
and then 

tm = tNe 
also 

Cm = G/Ne 


Substituting the above into the governing equation given by equation 


(16) yields . 


l 1 
Ne 5 (CmNe) + ae (UmUcGmNe) + ie (VmUcCmNe) = 


1 v2 8 eee 
TZ V¥m (GmNc) + LePo 0Xm (PoP mM) 
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simplifying 


P) U-N Pp) U.N 0 Ne d 
NE ee ES Sig mbm) +E ayes Wm) = TE YVR + ES 








which becomes 


oGm | nee O(umSm) ace O(Vm$m) Vv 9 dP'm 

Jin TLNe Okm *DeNo 8m Nez Vam tN, 136) 
but 

= —F = 1 

NZL. ~ NZ 


Thus equation (36) becomes 





atm , Uc _(aumbm) , avmém))_ __v__ Ue yor, 20m 
dtm Sides 0Xm  Oym ce aie Uc *m>™ ™ Oxm 
or 


Gm ir. os + Sob) Fy Deas OP'm 


otm n° dym |" Re Vmom+3y ea 


where Fy and Re are as defined above. 


The density diffusion equation is similarly nondimensional- 


ized as 


, , = Pe 
Ps one) yet ev Bo wp! +v SE 
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which becomes 


Ne > ot (PoP'm) + 7 - sc (UmUcPoP'm) + F- os (VmUcPoP'm = 





_UcVm _d v 2 es 
Co Sym PoPm) + TZVA (Popa) + Tear (PoP m) 


The above equation can be simplified to yield 


OP'm m | UcPo a 








; U d 
PONc St + Le Oxm UmP'n) + og Jym VmPm) = 
-U a0 Vv V0 07D 
= a. Vn = + AB Po Tz Vein + iz T 


Introducing Fy and Re as before and recalling that N2 = a one has 


OP'm O(Ump'm) , 2m’) OPm 
dtm * *v — Om * Oym <a Pym Oyen. 
Fy 0-0 m 

Re va P'mt+ 32 =| 


(38) 


Now as developed above in equation (34) 


7: N2_ NL 
2=- =—— 2= = 
n o¥m where n N Fi 





+= 
Substituting for = in equation (38) yields the nondimensionalized 
density diffusion equation 


ap’ a(ump'm) | 2(Vme'm) Fy op 
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Equations (37) and (39) are valid when Fy<<1l and buoyancy 


dominates the flow. When Fy approaches zero, the equations that 


result from equations (37) and (39) describe the propagation of linear 


internal waves. The buoyant scaled forms of the governing equations 


are of interest to the present investigation because, for submerged 


bodies of naval interest, Fy is typically about 0.001 and thus signifi- 


cantly less than one. 


Having established the buoyant scale as the form of interest, 


the boundary conditions can also be nondimensionalized using the 


same technique. From equation (21), 


d 
A udl ay aty =n 


which becomes 


Nm = N/Le 
Simplifying the above yields 


; onm Uc Nm Uc 


Otm * Neko U™Oxm ~ Nelo “™ 
Introducing Fy as before 


on on 
=. t Fy Um 5, = Fwm at Ym =m 
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(40) 


and 


Om = nae (nondimensional velocity potential) 
Cc 


Equation (25) can be nondimensionalized as follows: 


0D “ 
at a + en = 0 : 
becomes 


O(SmNcL2) U2q2 
Nese + 3 + Blot = 0 (41) 


where 
Pe? 2 
Gm = Um + Vn 


Further simplifying equation (41) yields 


Om ~ Be aes £m 


dtm ~ N22 2 * NLC 0 (42) 
but 

za - ¢ 

NZL = ] and N2LZ Jak 


Thus equation (42) becomes 


a oe 
st Fe 3 + 1m = 0 at Ym = "im 





32 





For the purposes of the present investigation, the flow will be 
assumed to be inviscid, i.e., the viscous diffusion will be ignored to a 
first order approximation. The nondimensionalized forms of the 
governing equations and boundary conditions are summarized in 


igure 1. 
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Tl. NUMERICAL SOLUTION TECHNIQUES 


A. INTRODUCTION 

The significant difference between this study and the work of 
previous investigators is the introduction of a deformable free surface 
and the resulting complex interplay between kinetic and potential 
energies. The computational domain thus has an unknown boundary 
on which a double condition must be imposed as previously outlined. 
This complexity, as well as several other specific features of this 
problem, directly influences the ability of a numerical method to con- 
verge to an adequate solution. 

The governing differential equations and boundary conditions are 
both nonlinear. Sarpkaya and Issacson [Ref. 23] note that for small 
amplitude waves (amplitude << wavelength), the boundary conditions 
at the free surface may be linearized. This approach, although possibly 
valid in the vicinity of surface striations, would not be valid for the scar 
front where observed light diffraction patterns [Ref. 17] indicate very 
small radii of curvature. Additionally, Haussling and Coleman [Ref. 24] 
have demonstrated numerically the importance of nonlinear terms in 
solving the free surface problem in the vicinity of the generation 
source. This is the case encountered in the region of the scar front. 

The rate of change of the free surface deformations also tends to 
be slow. The scar pattern develops as the vortex pair rises to its 


maximum height and then the scar is trapped by and slaved to the 
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vortex beneath the surface. The scars and striations are thus not small 
amplitude surface waves propagating independently from their 
generator, but are in fact local disturbances whose generation, growth, 
and propagation are controlled by the vortices. This aspect, observed 
experimentally by Gray [Ref. 17], separates this problem from other 
research in the field of deformable free surfaces. Previous investiga- 
tors, applying numerical methods to such surfaces, have dealt almost 
exclusively with surface waves generated by a moving body (surfaced or 
submerged) and propagating away independently. 

The numerical methods reviewed for solution of this problem can 
be broadly categorized into Finite Differences, Finite Elements, 
Boundary Integral Equations, and Hybrid Methods which employ com- 
binations of the others. For the purposes of this study, such catego- 
rization is based on the manner in which the governing equations are 
tackled in the computational domain. It is noted that all of the 
methods reviewed utilize some form of finite difference scheme with 
respect to time. The applicability of each method and the Poeneiated 
advantages and disadvantages are discussed and summarized at the 


end of this chapter. 


B. FINITE DIFFERENCE TECHNIQUES 

The use of Finite Difference techniques in solving partial differen- 
tial equations has been well established and successfully implemented 
by several investigators. For the case of a nondeforming free surface, 


the governing equations are solved in the computational domain using 
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of . a(ug) | atv )_ vv2t _ £ 9p" (30) 


and 


' om Ub Pa 


and the Biot-Savart integrals 


Hae va fs, ') dx'dy’ (31) 
A 
ae si ' ’ ' . 
vease | (x -x') ao dx'dy (32) 
x | 


where 


r2 = (x -x'))+y-y’ 
x', y': position of the vorticity 


x,y: point where u and v are to be calculated 


In nondimensional form, the buoyantly scaled forms of equations (31) 


and (32) become 


Sor 


(Yim — Ym) Gm ('m. Ym) dX mdym L,.Nec 
Hn bimym) = JOE eee 
A 


(Xi — Xp) Gm (Xm. Ym) dXmdym L,Ne 
Vm (Xm.Ym) = —— Sarl Uc 
A 


<< 


Upon simplification, the equations become 


im — Ym) (X'm. Ym) dxX'm dy’ 
Rye neal 7 _ Ym Se Xm dYm (33) 
A 
6% — Xin) Cm (X'm: Yin) OX'm dy'm “we 


Fy Vm (Xm.Ym) = | Dare 
A 


where parameters are nondimensionalized as before. 


The governing equations can then also be rewritten as 


Gm , oly Um Gm) O(Fy Vm Gm) _ Fy v2 OP'm 
dtm ~ OXm 0Ym ~ Re /m Gm + Xm 
and | 
pm | o(Fy Um Pm) | O(Fy Vm P'm) _ 9 Fy 0. 020m 
be Se ay ae -=FyVmn nce V Pm * oy2— 


where Fyum and Fyvpm are calculated directly using the nondimension- 


alized forms of the Biot-Savart Integrals equations (33) and (34). 
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The treatment of the governing differential equations in a stream 
function formulation is thus relatively straightforward. The complica- 
tion is associated with the determination of the free surface location in 
conjunction with the solution of the governing equations, all preferably 
simultaneously, as noted by Yeung [Ref. 25]. 

Finite Difference methods are most suitable, or at least simplest 
to implement, for a boundary geometry that is rectilinear. The 
deformable free surface requires the use of “irregular stars” for differ- 
ence schemes near the fluid boundary. These “irregular stars” are of 
an inferior accuracy when compared to the remainder of the grid. 
Thus, either the mesh must be refined or the difference formula 
changed to a higher order if accuracy is to be maintained near the free 
surface and especially in the vicinity of the scars. Yeung [Ref. 25] 
notes that, with the location of the free surface unknown a priori, we 
have the unfortunate situation that the regions that demand the great- 
est accuracy are precisely those where it is hardest to achieve. Also, 
since the shape of the free surface affects the migration of the vortex 
pair, a loss of accuracy in determining the free surface will be 
reflected in the determination of the location of the vortex pair. This 
complication is one not included in the investigations reviewed below 
where the generating body moves independently of the free surface 
location. 

Conformal transformations can simplify one aspect of the problem 
by transforming the real geometry with a deformed free surface to a — 


rectilinear mesh. The boundary conditions are thus simplified but the 
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resulting field equations are increased in complexity. Haussling and 
Coleman have successfully utilized transformation functions of this 
form to generate a time independent computational region in which 
nonlinear free surface boundary conditions are applied. 

The Finite Difference method primarily benefits from the direct 
solution. of the governing differential equations as opposed to- the 
development of intermediate integral relations required in the inte- 
gral boundary techniques. The cost of such an approach is in the 
computational intensity required to iteratively solve the finite differ- 
ence forms over the computational domain. Haussling and Coleman 
have demonstrated the use of a successive over relaxation technique to 
solve steady state nonlinear free surface boundary condition problems 
as discussed above. 

A modification of successive over relaxation, and one offering 
increased rates of convergence, has been demonstrated by Brandt, 
Dendy, and Ruppel [Ref. 26]. Their technique utilizes a multigrid 
solution which solves for low-frequency components of error on a 
course grid where the calculation is relatively inexpensive, and high- 
frequency components of error on a fine grid where successive over- 
relaxation is efficient. 

Theodussiou and Sousa [Ref. 27] also utilized a modified grid sys- 
tem to speed convergence by staggering their grid structure such that 
pressure was defined at the center of the discretized control domain 
and velocity components were defined at the center of the control 


domain faces. This arrangement has the convenient feature that the 
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velocity components are stored at just the points at which they are 
required for the calculation of their advective contributions. The 
pressure gradients can be represented by central differences without 
inducing non-physical oscillations in the pressure distribution. Note, 
however, that both of these multigrid techniques add to the complex- 
ity of establishing a mesh which moves with the deforming free-sur- 
face. Both sets of authors suggest, however, that this added complex- 
ity is outweighed by the savings in ease of convergence. 

Ohring [Ref. 28] and Ohring and Telste [Ref. 29] have also partially 
circumvented the computational intensity of an iterative technique by 
directly solving the finite difference equations resulting from the 
Laplace Equation. The technique employed utilizes a fourth-order 
solver to diagonally decompose the resulting coefficient matrix. Taylor 
series approximations are also suggested for application to the nonlin- 
ear free surface ere conditions. It is suspected, however, that 
the computational benefits derived from such a technique will be lost 
when the coefficient matrix becomes nonlinear as would result from 
the governing equations in this problem. 

For completeness, as noted previously, finite differencing in time 
is required for all the numerical methods reviewed. Sarpkaya and his 
students [Refs. 17-20] have successfully employed an upwind-differ- 
encing scheme in time and verified it with experimental results. 
Yeung notes that the free surface conditions, being first order in time, 
can be used to advance the solution of the elevation and velocity 


potential on the free boundary. However, the difference form utilized 
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can dramatically affect the stability and thus a modified Euler method 
is suggested since it is unconditionally stable. Haussling and Coleman 
successfully utilized this technique for advancing their solution in 


time. 


C. FINITE ELEMENT METHODS 

The Finite Element and Finite Difference methods enare! the 
common feature that both attempt to solve the governing differential 
equations directly, differing only in methodology. The Finite Element 
Method is one based on the method of weighted residuals. The usual 
procedure consists of first subdividing the domain of interest into a 
mesh of finite-sized subregions, within each of which the solution is 
represented by some convenient choice of trial functions, usually 
polynomials. The trial functions are determined by substituting them 
into the governing equations and requiring the integrated error or 
residual based on certain weighing functions to vanish. An integration 
by parts is normally preformed to reduce the interelement continuity 
requirements of the trial functions and to incorporate the 
nonhomogeneous boundary conditions. The weighing functions, also 
known as test functions, can be chosen in a variety of ways. A “weak 
formulation,” such as the Galerkin Method, makes the space of the 
test functions identical to that of the trial functions. In contrast, a 
“strong formulation” is one based on the existence of a variational 
principle where a functional is made stationary. Specifics of finite 
element techniques may be found in Zienkiewicz [Ref. 30] and Dhalt 
and Touzot [Ref. 31]. 
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The nonlinear aspects of this problem affect the Finite Element 
Method in much the same way as the Finite Difference Method. The 
coefficient matrix of a representation such as 

[K] [U] = [B] 
where 

[K] = eoefficient matrix ~ 

[U] = matrix of unknowns 

[B] = forcing function matrix 
is neither symmetric nor independent of [U], as in the case of a linear 
problem. An iterative technique is thus required to solve this problem 
at a given instance in time. As a result, much of the advantage of 
reduced storage and computation time inherent to finite element for- 
mulations is lost. The choice of iterative techniques does not differ 
from that available to Finite Difference Methods and, thus, once the 
computational domain is discretized, no significant difference in 
problem solution is involved. 

The discretization of the computational domain is an area where 
the Finite Element Method does have distinct advantages. The intro- 
duction of curvilinear or isoparametric elements of higher order 
allows one to cope with any arbitrary boundary geometry with little 
loss of accuracy. Yeung notes that, although this particular advantage 
is reduced when Finite Difference Methods are used with conformal 
transformations, Finite Element Methods still retain superiority in 
flexibility, particularly in the case where varying size and shape 


elements are introduced to overcome local irregularities. Such 
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irregularities are exactly the problem involved in the vicinity of the 
deformable free surface. Curvilinear elements could be used to fit the 
free surface with finer elements incorporated in the vicinity of the 
scar front. This could be accomplished with little loss of accuracy and 
less complexity than that which results from the use of “irregular 
stars” in the Finite Difference Method. Note, however, that, since the 
scar front moves with time, an adaptive mesh refinement technique 
would be required to appropriately place the finer elements in the 
vicinity of the scar. Sarpkaya and Hiriart [Ref. 32] successfully utilized 
varying size elements in the vicinity of the free surface in conjunction 
with their moving net computation. Although the free surface in their 
case differed from that involved in this problem, the basic concept of a 
“flexible mesh” remains unchanged. 

Admittedly, the basic problem of determining the location of the 
‘nee Rentice in conjunction with solving the governing equations 
remains one of trial and error, regardless of the type of discretization 
employed. Larock and Taylor [Ref. 33] adjusted the free surface loca- 
tion to achieve tangency to the surface velocities calculated based on 
an assumed free surface position. The pressure boundary condition 
was then used as a check of this resulting location. Larock and Taylor 
note, however, that such a technique will not work well where high 
curvature or Froude numbers (Fy) less than one are involved, as is the 
case in the present investigation. As an alternative, Sarpkaya and 
Hiriart adjusted their free surface location to satisfy both the pressure 


and velocity boundary conditions simultaneously. 
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The Finite Element Method, when formulated using the varia- 
tional techniques, also benefits from the fact that the boundary condi- 
tions may be included as an integral part of the functional rather than 
dealt with separately. Bai and Yeung [Ref. 34] and others have suc- 
cessfully employed variational techniques to directly incorporate 


boundary conditions in the solution of linear water wave problems. 


D. BOUNDARY INTEGRAL EQUATION METHODS 

The term “Boundary Integral Equation Method” can be applied to 
a large group of numerical techniques that include Green’s Function 
formulations, Spectral Methods, and a boundary element application of 
the Finite Element Method. In all cases, the solution approach differs 
dramatically from that of Finite Differences and Finite Elements in 
that the governing equations are not attacked directly. Instead, the 
. problem is solved by satisfying boundary conditions on a discretized 
boundary and thus reducing the spatial dimensions of the problem by 
one. Physical quantities such as wave height and fluid pressures are 
required and solved for only on the boundaries. Interior data, although 
available based on the boundary solution, is not specifically required. 
Boundary Integral Equation Methods thus have the distinct advantage 
that only the physical quantities specifically desired on the free sur- 
face are required for the solution. 

Basically, there are two approaches to boundary integral formula- 
tions. Either an approximating function is chosen on the boundary 


that satisfies the governing equations in the domain and approximates 


4a 4 


the boundary conditions, or vice versa. In both cases, the solution 
technique can be further divided into indirect and direct methods. 

In indirect approaches, the fundamental solution is approximated 
on the boundary by a function with unknown coefficients. These coef- 
ficients are found by satisfying the boundary conditions. Distributed 
singularity methods, such as simple-source distributions, are an exam- 
ple of this approach. In this case, the fundamental solution is repre- 
sented on a discretized boundary by distributed simple sources of 
unknown strength at known locations. The strength of each source is 
then determined based on the solution of the boundary conditions. 

Direct methods find the fundamental solution through the use of 
Green’s Function formulations, which directly incorporate the 
governing equations. The resulting Green’s Function integrals are 
typically solved using quadrature or point kernel techniques. The 
main disadvantage in this approach is that, for the governing equations 
involved in this problem, there is no straightforward Green’s Function 
formulation that leads to fundamental solution. 

The main advantage of the Boundary Integral Equation Method is, 
then, the ability to directly discretize the free surface without a loss of 
accuracy. Information is thus obtained exactly where required. How- 
ever, Brebbia [Ref. 35] notes that this is not without the sacrifice of a 
symmetric coefficient matrix such as that which is common to the 
Finite Element Method. Yeung notes that, in general, the great 
reduction in the size of the matrix outweighs the added complexity of 


solving a nonsymmetric system of equations. 
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The disadvantages of Boundary Integral Equation Methods are 
directly associated with errors that are the result of discretization. 
This discretization plays such a large role in the formulation of an 
efficient boundary integral equation that Brebbia refers to these tech- 
niques as “Boundary Element Methods.” In all formulation tech- 
niques, -elements of one type or another are formed over whieh a 
fundamental solution must be approximated either by distributed sin- 
gularities, Green’s Function integrals, or trial functions for finite 
elements. 

Discretization errors result both from collocation errors and geo- 
metric surface errors. Collocation errors are primarily the result of 
leakage, as noted by Hunt [Ref. 36]. Since the boundary conditions are 
satisfied exactly only at discrete points, the remainder of the free sur- 
face is “porous” by comparison. Leakage errors are reduced by 
increasing the number of free surface elements and thus the number 
of discretization points. This, of course, is done at the expense of 
computational intensity. Geometric surface errors — from the 
approximation of a curved free surface by linear elements. This prob- 
lem becomes particularly noticeable in areas of high curvature, such as 
scar fronts. Higgins and Cokelet [Ref. 37] noted that this problem can 
be partially circumvented by using a Lagrangian description of 
“marked” particles on the free surface. These particles tended to 
concentrate in the regions of highest curvature, thus giving improved 


accuracy exactly where needed. 
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Finite Element Method formulations using boundary ieeiaaenis 
only can also circumvent some of the errors associated with the collo- 
cation and geometric surface by using curvilinear or higher order 
parametric one-dimensional elements. The problem then is one of 
finding suitable fundamental solutions that satisfy both the governing 
equations and continuity requirements between elements. - 

Finally, Boundary Integral Equation Methods are complicated by 
the necessity to include nonlinear terms in the boundary condition for 
the free surface. As previously noted, this is necessary in order to 
maintain accuracy in the vicinity of the scar front where very small 
radii of curvature occur. Several researchers have successfully applied 
various forms of the Boundary Integral Equation Method to Laplace | 
equations with linearized free surface conditions, but only a few have 
incorporated a nonlinear condition. Faltinsen [Ref. 38] incorporated 
nonlinear conditions by using the properties of the fluid particles on 
the free-surface at one instance in time to establish a new free surface 


location and step the solution forward in time. 


E. HYBRID METHODS 

The use of a different technique in different portions of the com- 
putational domain results in a hybrid numerical formulation. This type 
of formulation attempts to maximize the benefits of any one particular 
method by employing it only in regions where its accuracy remains 
high. A secondary objective is to reduce the computational intensity of 
the overall routine by using a coarser, less time-intensive technique in 


areas where accuracy is not of particular concern. This is exactly the 
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case in the regions far removed from the free surface. The usual 
approach is to take advantage of the availability of analytical solutions 
in regions where the flow geometry is relatively simple. This 
approach allows the number of mesh points to be reduced with a 
resulting decrease in required storage and computational intensity. A 
further advantage is that the analytical solutions can be chosen to 
permit a simple solution of radiation type boundary conditions, as 
noted by Yeung [Ref. 25]. 

The disadvantages of Hybrid Methods are associated with the 
necessity to properly match the solutions of different subdomains in 
the overall computational domain. This matching may result in 
numerical perturbations to the solution technique if not properly 
employed. Additionally, the advantage of reduced total grid points is 
somewhat offset by the added computations required to match 
solutions at the common boundaries. 

Yeung notes that “the more successful Hybrid methods have so far 
been restricted to linearized problems where analytical solutions in 
the exterior regions could be obtained without too much difficulty. In 
particular, treatment of steady flows in a uniform stream or time-har- 
monic flows with linearized boundary conditions have been quite well 
established.” [Ref. 25] 

Bai and Yeung [Ref. 34] utilized a finite element grid in the vicinity 
of the free surface and an eigenfunction solution in the exterior region 
to reduce the computational intensity of their routine. Chang and Pien 


[Ref. 39] used a source distribution to formulate a Boundary Integral 
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Equation Method near the free surface and a finite difference routine 
in the exterior regions to calculate the hydrodynamic forces on a body 
moving beneath a free surface. It should be noted, however, that both 
sets of investigators used a Laplace formulation with linearized bound- 
ary conditions to obtain their results. 

For.completeness, it is noted that Hybrid methods could be-con- 
sidered to include various numerical techniques that employ the same 
method throughout the computational domain but in different man- 
ners in each subdomain. Such is the case, for example, when a finite 
element method is employed with a varying grid and/or element type 


in various regions of the computational domain. 


F. NUMERICAL METHODS, CONCLUSIONS, AND SELECTION 

The synopsis of advantages and disadvantages for each numerical 
method listed in Table 1 provides a good source of information for 
making a wise choice of formulation to be used in this problem. The 
“best” method is one which will meet the goals of the investigation 
while maximizing the advantages in areas of particular concern. 

In this case, accuracy in the vicinity of the free surface is of par- 
ticular importance since ultimately it is the free surface shape which 
is desired. This accuracy must be obtained with full consideration of 
the necessity to include nonlinear boundary conditions while mini- 
mizing the computational intensity of an iterative process. Tuck, ina 
paper by Bai and Yeung, observed that “to a certain extent the ‘best’ 
method will always be that which appeals most to the person pro- 


gramming it, and hence that which gives him the greatest chance of 
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writing a successful program irrespective of efficiency. There is no 


less efficient program than one which does not work at all.” [Ref. 34] 


Given that a successful program can be developed (a pretty big given), 


it is then necessary to maximize both efficiency and accuracy. 


With full consideration of the concerns given above, the Boundary 


Integral- Equation Method surfaces as the method of choice for the 


following reasons: 


ie 


Zi 


Computational intensity is minimized by reducing the spatial 
dimensions by one. This aspect is particularly noticeable when 
consideration is given to the iterative requirements of this 
problem. This substantial improvement directly incorporates 
any similar advantage that can be gained by using a varying Finite 
Element Method or Hybrid Method. 


Accuracy is maintained at the free surface by incorporating the 
best aspects of Finite Element Methods and a Lagrangian 
description of marked particles. This direct discretization of 
the free surface provides both accuracy where needed most and 
an exact solution of boundary conditions at the marked points. 


Storage requirements are drastically reduced since only physical 
quantities at the free surface are required or needed. 


The implementation of the Boundary Integral Equation Method 


through the use of distributed vortices will be outlined in detail in the 


following section. 
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IV. DISTRIBUTED VORTEX MODEL 


A. INTRODUCTION 

Discrete vortices, with or without a core, or vortex sheets have 
been used as boundary elements to simulate separated and unsepa- 
rated flows. The method consists of the determination of the appro- 
priate strength of the vortices at each time interval through the use of 
the governing equations. The vortices are then convected to their 
new positions and the process is repeated. In spite of its simplicity, 
the distributed vortex model presents several difficulties, all of which 
are related to discretization and the use of vortices. The evaluation of 
the governing integral equation cannot accurately be accomplished 
merely by applying a standard integration formula. A vortex sheet or a 
string of point vortices is unstable to small sinusoidal disturbances of 
any wavelength. This phenomenon of Helmholtz instability persists in 
curved nonuniform vortex sheets, at least for short waves, unless the 
sheet is rapidly stretching. In other words, the round-off and trunca- 
tion errors are rapidly amplified to cause the chaotic motion which 
often ruins practical calculation. If it is granted that it is the growth of 
short waves which can ruin calculations with vortex sheets, it is sensi- 
ble to consider ways of removing the instability. This is because the 
instability is introduced by the step of replacing a shear layer of small, 
but finite, thickness by a vortex sheet. One could give up the vortex 


sheet approximation and return to the computation of the evolution of 
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a thin layer. This is without doubt the most satisfactory procedure, but 
it involves much more computation. 

An alternative approach is to modify the governing differential 
equation to allow for finite thickness but the resulting equation, while 
only a little more complicated than the original equation, has not 
proved amenable to computation. = 

Another possibility is to apply a linear smoothing formula, such as 
that introduced by Longuet-Higgins and Cokelet [Ref. 37] in their 
work on nonlinear water waves. The subsequent applications, includ- 
ing the one discussed herein, have shown that chaotic motion sets in 
sooner or later regardless of the smoothing. The repositioning tech- 
nique introduced by Fink and Soh [Ref. 40] and used subsequently by 
Sarpkaya and Shoaff [Ref. 41] removes the most unstable mode and 
reduces the growth of the higher modes of Helmholtz instability. 
However, it does not prevent the growth of spurious waves along the 
vortex sheet. The disadvantage of the smoothing, either through the 
use of a numerical filter or through the repositioning of the vortices, is 
that it is not clear in general what the relationship is between the 
results achieved and the unknown exact solution. The problem may- 
not possess a Solution for all time, and in this case the use of smooth- 
ing could yield an acceptable-looking solution where none in fact 
exists. Alternatively, the solution arrived at through smoothing may 


not be even close to the exact solution. 
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B. APPLICATION OF THE DISTRIBUTED VORTEX METHOD 

The liquid surface was represented by a number of point vortices. 
Their positions were symmetrical with respect to the x axis, situated 
on the free surface. However, the strength of a vortex at (x, y, nor- 
malized by bo) was opposite to that of a vortex (image vortex) at (-x, y). 
From a mathematical point of view, one would like to have the vortices 
extend from -< to +e and the two trailing vortices originate at (1/2, 
-coo) and at -1/2, --). This is impossible from a numerical point of 
view. Observations have shown that the free surface rises in a rela- 
tively small region directly above the trailing vortices. The remainder 
of the free surface remains undisturbed. These observations and 
several sample calculations led to the conclusion that the free surface 
can be restricted to a region extending from x = -10 to x = +10. 
Furthermore, the trailing vortices are assumed to originate at a depth 
of y = -5. 

The complex function representing the two trailing vortices 


below the free surface and the vortices on the free surface is given by 


_ We To Pm 
ii a In (Z - Zo) + BF ln (z + Zale In (Z - 2) 


ilm : 
= a (Z + Zi) 
in which the strengths of the free-surface vortices are normalized by 


the strength of the trailing vortex. All coordinates are normalized by 


by. Then the velocities, normalized by Vp = (T9o/2nbg) anywhere in the 


fluid medium is given by 
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meal 2M (z-Zo) 2% (z+Z_) 2% (2-21) 2% (747%)) 











(35) 
The normalized boundary condition at the free surface is given by 


Di, 2 + FO -) 


Thus, at least theoretically speaking, one can calculate the poten- 
tial function @ from the complex velocity potential w (the real part of 
w), the velocity of the vortices from equation (35), the elevation of the 
free surface from equation (36) for a given Froude number Fv (Fv = 
Vo/V bo), and trace the evolution of the free surface as a function of 
time. This relatively simple-sounding procedure is anything but sim- 
ple, primarily because of the fact that the numerical instabilities even- 
tually lead to large-scale instabilities, as noted in the introduction to 
this section. The use of various smoothing techniques was a viable 
option and, in fact, was tried at various stages of the calculations. The 
results have shown that the growth of the instabilities increases with 
decreasing Froude number. Neither the use of smoothing techniques 
(e.g., the Longuet-Higgins technique) nor the use of vortex sheets, 
vice discrete vortices, can eliminate the eventual development of a 
chaotic behavior on the free surface. It is because of this reason that 
the use of smoothing was ruled out and the calculations were 
restricted to a relatively high Froude number (Fv = 1.125). 

The specific details of the numerical calculations are as follows: 
(1) assign the position of the vortices; (2) find the strength of the vor- 


tices and the velocity potential assuming the free surface to be rigid; 


of 


(3) advance the position of the trailing vortices for a time interval t 
(0.01 in the calculations); (4) recalculate the velocity of the surface 
vortices and advance them forward in time for a time interval t; (5) 
calculate the strength of the free surface vortices by iteration until the 
free surface condition is satisfied; (6) calculate the velocity of the sur- 
face and trailing vortices and advance their positions for a time inter- 
val t; and (7) repeat the calculations by returning to step (5). 

The procedure described above is relatively simple and does not 
require excessive computer time (about 20 minutes on IBM 3033). 
However, the free surface eventually does become chaotic. Some of 
the instabilities can be alleviated without smoothing through a judi- 
cious selection of the initial position of the surface vortices. Numerous 
calculations have shown that the vortices near the y axis begin to slide 
sideways as the free surface (or the vortex sheet) stretches. Conse- 
quently, the thinly populated regions of the sheet do not yield a 
smooth free surface and the flow begins to leak between the vortices. 
Also, the regions where the free surface is depressed (where the scars 
develop) become overpopulated with vortices, leading to the growth of 
short-wavelength Helmholtz instability. This problem can easily be 
alleviated by packing the vortices more densely near the y axis at the 
Start of the calculations so that the entire surface becomes more or 
less uniformly represented at later times. This simple techn has 
been used in the results presented herein. An exponential function 


was employed to assign the initial vortex spacings so that near the 
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y axis the vortices were separated by a distance of about 0.03 and by a 


distance of about 0.25 towards the end of the vortex sheet. 


C. RESULTS OF THE NUMERICAL CALCULATIONS 

The evolution of the free surface with time is shown in Figures 2.1 
through 2.27 for Fv = 1.125. The normalized. time in these plots is 
given by T = (Vo t/bo ~ do/bo) where do is the initial depth of the 
trailing vortices. The time T = O corresponds to that at which the 
trailing vortices would have reached the undisturbed free surface had 
they continued to move with their initial mutual induction velocity. 
The use of other normalized times is not suitable since they tend to 
depend on the initial position of the trailing vortices. mi 

Figures 2.1 through 2.27 show that the free surface directly above 
the vortices rises rapidly while the adjacent portions of the surface 
depress and form two strong scars. Unlike the rigid surface case, 
where the trailing vortices continue to move towards right and left at a 
depth of about y = -0.5, the trailing vortices below the deforming sur- 
face almost come to rest at a point slightly above the free surface. 
Furthermore, their initial spacing remains nearly constant. This sug- 
gests that the Kelvin oval formed by the trailing vortices pushes the 
free surface up as if it were rigid during its upward migration. Evi- 
dently, this finding is based on the assumption that the trailing vor- 
tices are point vortices and are not subjected to viscous and turbulent 
diffusion. In reality, the vortices quickly become turbulent, their core 
radius increases, and the vorticity is diffused over an ever-increasing 


area with the passage of time. The amount of diffusion depends on the 


og 


initial position of the trailing vortices. Consequently, the trailing vor- 
tices emanating from large depths diffuse over a larger area relative to 
those emanating from smaller depths. There is at present no suitable 
mathematical means to deal with the turbulent diffusion of such vor- 
tices. The best one can do is to generate the trailing vortices at 
depths sufficiently close to the free surface to prevent excessive dif- 
fusion and yet sufficiently far so that the free surface remains 
undeformed at the start of the motion. Extensive calculations and 
experimental observations have shown that the free surface does not 
begin to deform until the trailing vortices reach a depth of about y = 
-1. Thus, it is perfectly safe to place the trailing vortices at y = -3 at 


— 


the start of the calculations. 
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V. EXPERIMENTS 


A. EXPERIMENTAL APPARATUS AND PROCEDURES 

Experiments were conducted in a water basin. It consisted of a 
12-foot-long, 3-foot-wide tank with aluminum bottom and plexiglass 
walls (see Figure 3). Additional equipment consisted of plumbing for 
the filling and emptying of the tank, a collimated light source, and the 
dye system for flow visualization. | 

The two-dimensional trailing vortices were generated through 
the use of two counter-rotating plates (see Figure 3). The mechanism 
rotating the plates was located at the bottom of the tank and ‘below a 
plexiglass plate spanning the entire tank. In other words, the motion 
of the driving system did not disturb the flow above the plexiglass. 
The mechanism was actuated by releasing a weight attached to the 
common driving shaft. Fluorescent dye was introduced slowly into the 
region between the plates through the use of two holes connected to 
two dye reservoirs. The reason for the use of two reservoirs was that 


different colors of dye can be used to visualize each trailing vortex. 


B. PROCEDURES 

Experiments were initiated by filling the tank to a suitable level, 
removing any air bubbles, bringing the plates to their full open posi- 
tions, waiting for a sufficient period of time for the fluid to come to 


rest, introducing the neutrally buoyant dye, actuating the light source 
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and the video system, and initiating the rotation of the plates by sud- 
denly releasing the load attached to the driving shaft. The plates 
rotated smoothly and then came to rest on the plexiglass. In other 
words, the plates generating the vortices literally “disappeared.” 
Thus, no other vortices were generated. It is a well-known fact that 
this is not the case with piston-generated vortices. When the piston 
stops, two additional vortices (from a two-dimensional piston) or 
another vortex ring (from an axisymmetric piston) are generated. The 
mechanism used in the present investigation effectively prevents the 
generation of secondary vortices by simply disappearing. 

The trailing vortices rise under their mutual induction velocity, 
quickly give rise to a Kelvin oval, and smoothly approach the free sur- 
face. When the vortices are at a distance of about 1.0 from the free 
surface, the free surface begins to rise and forms a smooth hump, bor- 
dered by two scars. For the Froude numbers achievable (maximum 
0.35), the trailing vortices turn quickly to right and left (as if they 
were approaching a rigid surface), iit scar front moves in the respec- 
tive directions ahead of the vortex, and the hump between the vortices 
begins to recede. During the later stages of the motion, the scars are 
slaved to the vortices and continue to move ahead and in the direction 
of the vortices. 

af 4.1 through 4.12 show the evolution of the free surface at 
various times T for a Froude number of Fv = 0.35. The vortex centers 
are clearly visible in these pictures because of the additional assistance 


provided by the dissolved air bubbles to the flow visualization efforts. 
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VI. DISCUSSION OF RESULTS 


Maen coliputer code based on the boundary element methods 
through the use of the discrete vortices has been used to carry out 
calculations for various Froude numbers in order to predict the evolu- 
tion of the free surface deformation. The type of the vortex distribu- 
tion, time increment, and the iteration techniques have been varied 
within limits to explore the differences in the motion of the free sur- 
face. The results have shown that the calculations are fairly stable at 
relatively high Froude numbers. However, at relatively small Froude 
numbers, the Kelvin-Helmholtz instability sets in quickly and the free 
surface exhibits highly irregular shapes. It appears that either the 
integration procedures have to be refined or more sophisticated vor- 
tex sheets have to be used to calculate the small deformation of the 
free surface at small Froude numbers. 

The results obtained with a Froude number of Fv = 1.125 are 
shown in Figures 2.1 through 2.27. The most striking feature of these 
results is that the free surface rises rapidly to a height of about 1.25 
above the mean water level and captures the vortex pair. It has not 
been possible to carry out the calculations to larger times. ‘The fact 
Sould be kept in mind that the point vortices used in the numerical 
calculations may have very little resemblance to the real vortices by 


the time they rise to the mean water level. In fact, the experiments 
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show that the vortices become quickly turbulent. There is at present 
no possibility of incorporating into a numerical analysis the motion of a 
turbulent vortex. 

Experiments were carried out at a maximum Froude number oF 
0.35. The evolution of the free surface is shown in Figures 4.1 through 
4.12. The shape of the free surface at the time of maximum rise is 
plotted in Figure 5. Figures 4 and 5 show that the free surface rises to 
a maximum height of about 0.25, develops two strong scars, and then 
subsides as the vortices begin to move parallel to the free surface. The 
scar front is slightly ahead of the vortex center and reduces to a small 
depression at larger times. The results presented in Figures 4 and 5 
should form the basis of comparison for future numerical efforts at the 
corresponding Froude numbers. Additional calculations are underway 
with more sophisticated vortex sheets. The results presented herem 
are extremely encouraging and are expected to lead to a better under- 


standing of this extremely complex and challenging phenomenon. 
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VII. CONCLUSIONS 


The investigation reported herein warranted the following 


conclusions: 


1. 


The basic equations governing the motion of vortices in homo- 
geneous and density-stratified media have been developed and 
expressed in terms of dynamic and buoyant scaling laws; 


A novel experimental technique has been devised to generate a 
nearly two-dimensional vortex pair rising toward a free surface; 


Migration of a vortex pair toward a free surface gives rise to a 
bulge and two scars; 


The maximum rise at the free surface occurs when the vortex 
pair is at a depth of about 0.5; 


In the range of Froude numbers encountered in the experi- 
ments, the rise of the free surface is always followed by a fall as 
the vortices begin to move parallel to the free surface at large 
times; 


The numerical analysis of the free surface deformation through 
the use of discrete vortices leads to instabilities at low Froude 
numbers. These instabilities could have been removed with a 
suitable numerical filter or smoothing technique. However, the 
use of a relatively subjective smoothing procedure has been ruled 
out. The calculations at a Froude number of 1.125 have yielded 
fairly stable solutions and provided the basis for comparison with 
future experiments. 
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Figure 2.14 Velocity Field for T* = -1.10 
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Figure 2.15 Velocity Field for T* = -1.00 
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Figure 2.16 Velocity Field for T* 
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Figure 2.19 Velocity Field for T* = -0.60 
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Figure 2.22 Velocity Field for T* = -0.30 
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Figure 2.26 Velocity Field for T* = 0.10 
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Figure 2.27 Velocity Field for T* = 0.20 
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gure 4.6 Vortex Motion at T* 
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Vortex Motion at T* = -0.09 


Figure 4.8 
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Figure 4.10 Vortex Motion at T* 
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